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Abstract 

We discuss the properties of an analytical solution for waves in radiating fluids, 
with a view towards its implementation as a quantitative test of radiation hydrody- 
namics codes. A homogeneous radiating fluid in local thermodynamic equilibrium is 
periodically driven at the boundary of a one-dimensional domain, and the solution 
describes the propagation of the waves thus excited. Two modes are excited for a 
given driving frequency, generally referred to as a radiative acoustic wave and a 
radiative diffusion wave. While the analytical solution is well known, several fea- 
tures are highlighted here that require care during its numerical implementation. 
We compare the solution in a wide range of parameter space to a numerical in- 
tegration with a Lagrangian radiation hydrodynamics code. Our most significant 
observation is that flux-limited diffusion does not preserve causality for waves on a 
homogeneous background. 
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1 Introduction 

Analytical solutions for radiation hydrodynamics are difficult to obtain due to 
the complexity of the equations but are a powerful tool for testing complex 
multi-physics codes. One of the most useful simplifying assumptions for any 
set of nonlinear differential equations is to set up an equilibrium state and an- 
alyze small departures from that equilibrium. Since this approach retains most 
of the terms in the equations, it not only provides valuable physical insight 
but also serves as a comprehensive and sensitive test of numerical algorithms. 
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Many perturbation studies of radiation hydrodynamics have been performed; 
we follow closely Mihalas & Mihalas [Tf2] and Bogdan et al. [3], and refer 
the reader there for additional references. Despite the straightforward appli- 
cation of perturbation theory to the equations of radiation hydrodynamics, 
however, there appear to be few numerical tests of this type of solution in 
the literature (reference [3] is one example). Our goal here is to conduct a 
systematic comparison of such a solution with a numerical algorithm in a 
wide range of parameter space. The code that we use for comparison is the 
Lagrangian radiation hydrodynamics code Kull [3] . We begin in §|5] with an 
overview of our assumptions and the form the equations of radiation hydro- 
dynamics take under these assumptions. The failure of flux-limited diffusion 
to capture free-streaming radiation waves is highlighted in £j3l We discuss the 
analytical solution in £jU and compare our results with previous work in £j5j 
Numerical results are given in £0 and we summarize in §7J 



2 Assumptions and Equations 

We investigate perturbations from an equilibrium state of constant density and 
temperature with zero velocity and zero radiation flux. In addition to drop- 
ping terms that are higher than linear order in the perturbation amplitude, 
we further simplify the equations by making the following standard assump- 
tions: 1) the material fluid is an ideal gas, 2} the material and radiation are 
in local thermodynamic equilibrium (LTEj^J, 3) the opacity is independent 
of frequency, and 4) scattering is negligible. The common assumption of an 
opacity that is also independent of temperature and density is not strictly 
necessary for a perturbation analysis; one can easily show that variations in 
the opacity due to density and temperature perturbations give rise to terms 
that are higher than linear order in the analysis. 

The above assumptions must be supplemented with a prescription for the 
configuration of the radiation field. One approach is to solve the radiation 
transport equation directly, making some assumption for the angular distri- 
bution of the radiation [3]. Alternatively, one can calculate angular moments 
of the transport equation and invoke a prescription for closing the moment 
equations. A commonly employed closure scheme is the Eddington approx- 
imation, which assumes that the radiation stress is isotropic and given by 
P = (E/3)I, where E is the radiation energy density and I is the unit tensor. 
This is the approach taken by, for example, Mihalas & Mihalas [Tl[2] . 



Note that a common temperature for the material and radiation only applies 
to the equilibrium state; the material and radiation temperature perturbations are 
allowed to differ. 
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Both of these approaches are numerically expensive, however, due to the large 
dynamic range between the length and time scales of the material and radi- 
ation. The disparity in time scales can be alleviated somewhat by invoking 
the diffusion approximation, which assumes that the time dependence of the 
radiation flux is negligible. Since this can result in a superluminal flux of radi- 
ation energy, numerical calculations typically employ some type of flux-limited 
diffusion, which gives one the computational advantages of the diffusion ap- 
proximation while preventing the flux from becoming unphysical. 

As we discuss in the following section, however, all flux limiters reduce to the 
diffusion limit for linear perturbations. As a result, our numerical calculations 
are in the diffusion limit, although we discuss the analytical solution under the 
Eddington approximation for comparison with previous work. The equations 
of radiation hydrodynamics under the Eddington approximation in a frame 
comoving with the fluicO] are 



Dp 
Dt 



-pV ■ v, 



P^ = -Vp+*F, (2) 
Dt c 

= - 7 pV ■ v + c X (l -1)(E- a B T A ) , (3) 

^ = _ V • F - ~EV ■ v + c X (a B T A - E) , (4) 

"TiT = ~l VE - ~ FV - V -X F > ( 5 ) 

c Dt 3 c 



where p, p and T are the material density, pressure and temperature, re- 
spectively, v is the fluid velocity, F is the radiation flux, is the radiation 
constant, c is the speed of light, and x is the absorption opacity in units of 
inverse length. 



The perturbed form of the above equations is 

V • Sv = 0, (6) 



d_ 5p 
dt Ipo, 



7(9 (5v\ ( 5p\ 167r% SF 

adt\7 + \p~ ~4T r T)4E = ' U 



2 Since our perturbation analysis implies a constant opacity, the results are inde- 
pendent of the choice of reference frame. 
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dt \po Po) \T T J K J 

d (ST r \ „ / SF 1 \ (ST 5T r \ 



o / 

where T r is the radiation temperature defined via E = asT^ (E = obTq), 
a = (7Po/po) 1//2 is the material sound speed, and the dimensionless ratio 

r _ ( 7 - l)a B T* 

governs the coupling between the radiation and the material; it is proportional 
to the ratio of their energy densities. The subscript zero denotes an equilibrium 
quantity. For an ideal gas, the material pressure perturbation is given by 

s JL = s .L + S -l. (12) 

Po T po 

Under the diffusion approximation, the perturbed radiation energy and mo- 
mentum equations Qj and ffTUj) reduce to 

d (8T r \ c„ 2 (ST r \ 1 r (ST ST r \ . . 



and 



3 Breakdown of Flux-Limited Diffusion 



The radiation flux under the flux-limited diffusion approximation has the form 
(in the comoving frame) 

F = -—VE, (15) 

X 
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where A is the flux limiter, designed to reduce to 1/3 in optically thick regions 
and xE/\^E\ in optically thin regions. Flux limiters can take various forms 
but are generally nonlinear functions of 

R - (16) 



The diffusion limit corresponds to R <C 1 . For an inhomogeneous medium such 
as an atmosphere, |"Vi^| ~ E/L, where L is the characteristic length scale of 
the medium. The magnitude of R then depends only upon the optical depth 
r L = xL. In the case of small perturbations on a homogeneous background, 
however, 

„ k5E 

(17) 



where k is the wave number of the perturbation, so that R depends upon both 
a perturbation optical depth r k = x/k and a perturbation amplitude. A linear 
solution clearly requires R <C 1, even when r^T 1 ^> 1. All flux limiters there- 
fore operate in the diffusion limit for linear perturbations on a homogeneous 
background with zero mean flux. In regions of the flow where the diffusion 
approximation breaks down, this results in superluminal propagation speeds. 

This behavior, while somewhat counter-intuitive, is to be expected since flux- 
limited diffusion is not designed to follow wave fronts. The flux for large R 
(the free-streaming limit) is designed to reduce to F = cEn, where n is a 
unit vector in the direction of propagation; this constant flux can be viewed 
as a phase-averaged wave amplitude multiplied by a group velocity (i.e., a 
phase-averaged Poynting flux). Capturing the wave oscillations themselves 
clearly requires retaining the relevant time dependent and gradient terms in 
the equations; for free-streaming radiation, these are precisely the terms that 
are neglected in flux-limited diffusion (see, e.g., Levermore & Pomraning [B] 
equation [14]). 



4 Analytical Solution 



Perturbations on a homogeneous background are naturally decomposed in 
terms of Fourier modes, whose space-time dependence is exp(iu)t — ik ■ x). 
As discussed in Bogdan et al. [3], these plane wave solutions to the pertur- 
bation equations can be decoupled into modes parallel and perpendicular to 
the velocity. The transverse modes are akin to viscous shear modes in hy- 
drodynamics and we will not discuss them further here. For a system driven 
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at a constant frequency u, the longitudinal modes give rise to a fourth-order 
dispersion relation for the wave number: 

c 4 r fc - 4 + c 2 r,7 2 + Co = 0, (18) 



where the coefficients take different forms depending upon the approximation 
being used. For details on the derivation of the results in this section, see 
Appendix |A] 



4-1 Eddington Approximation 



The coefficients of the dispersion relation (fT8|) under the Eddington approxi- 
mation are 

c 4 = 1 — i!6rr c , (19) 



c 2 = 3(1 + it' 1 ) 2 - r-\l - zl6 7 rr c 

/ 1 IQ'-yr + it 

+ 16r 5 + Sir' 1 + 



-1 ' 

c 



3[7-i] y 

c = -3r- 2 (l + it- 1 + 16 7 r) ( 1 + it' 1 



16ra 2 
3[7 - l]c 



2 ' 



(20) 
(21) 



where 



Ta = ^ (22) 



and 



(23) 



As discussed in Bogdan et al. [3], the solutions to the dispersion relation 
(Tl8l) have a simple form in most of parameter space. In Figures [1] and [2] we 
reproduce Figure 4 of Bogdan et al. [3] under the Eddington approximation;. 
The leading order solution in the regions defined by Figured] for the radiative 
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acoustic mode is 
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These all have the form u/k a = v p (l + ie), where the phase velocity v p is 
a in regions a and b, the isothermal sound speed a/^/7 in regions c-e, the 

radiative sound speed J^P/ (3p) in region f (where P = E/3 is the radiation 

pressure), and c/ a/3 in region g. Except near the region borders, the radiative 
acoustic wave is weakly damped (its damping length is much greater than its 
wavelength) . 

The leading order solution for the radiative diffusion mode (Figure W) is 
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(25) 



This wave is strongly damped everywhere except region H and the r c < 1 
portion of region G. In regions A-D the damping length is on the order of the 
perturbation wavelength, and in regions E and F and the r c ^> 1 portion of 
region G it is much greater than a wavelength. The phase speeds are tJ p Ca 
in regions A and D and a portion of region E, a < v p < c in regions B, C and 
F and a portion of region E and v p = cj V3 in regions G and H. 
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4-2 Diffusion Approximation 

Under the diffusion approximation, the solution for the radiative acoustic wave 
remains the same everywhere except region g, where it can be seen from expres- 



the final term in cq becomes important. As noted by Bogdan et al. [3], how- 
ever, a consistent treatment of this region of parameter space would require 
relativistic physics, since the radiation energy density is greater than the rest 
mass energy density of the material. 

The solution for the radiative diffusion wave remains the same everywhere 
except regions G and H. Region C extends into region H and the portion 
of region G for which r c <C 1, so that the mode remains diffusive rather 
than becoming free-streaming. Its phase speed is super luminal and increases 
without bound as the driving frequency is increased!^ Its damping length is 
on the order of a wave length. In the portion of region G for which r c ^> 1, 
the solution is given by 



so that the mode is free-streaming with a phase speed (and group speed) ~ 1.2c 
and a damping length half as long as under the Eddington approximation. 



5 Comparison with Previous Work 

5.1 Mihalas & Mihalas 

Our coefficients (TT9T) - (T2T]) are somewhat different from those given by Mihalas 

6 Mihalas [TH2] due to their neglect of the velocity dependent term in the 
radiation momentum equation f[TUl)f^1 Comparison with expression (3.12) of 
Mihalas & Mihalas pQ reveals three differences: the term 3ir c _1 in c<i and the 
final term oc ra 2 /c 2 in Co are missing from the dispersion relation of Mihalas & 
Mihalas pQ, and they have an additional factor of l+ir~ l multiplying the final 

3 This is also true for the dispersion relation of Mihalas & Mihalas |l|2j ; see $5) 

4 Since diffusive modes have u> oc k 2 , their group speed is twice their phase speed. 

5 This term arises due to the Doppler shift between the comoving and laboratory 
frames. Mihalas & Mihalas [2] refer to it as an acceleration term since in the co- 
moving frame it appears as a time derivative of the velocity. 





(26) 
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terms in parentheses in C2 which cancels out in a self-consistent treatment. The 
only one of these differences that appears to be significant is the final term of 
c ; as we discussed in §4.21 this term is essential for limiting the phase speed 
to less than the speed of light at sufficiently large values of r. 



5.2 Bogdan, Knoelker, MacGregor & Kim 

Bogdan et al. [3] attempt to capture both the optically thick and optically 
thin regimes by explicitly calculating angular moments of the perturbed in- 
tensity, obtained directly from the perturbed transport equation. They obtain 
a transcendental equation for the wave number that reduces to a quadratic 
equation in both the optically thick and optically thin regimes. Remarkably, 
despite significant differences between the coefficients of their dispersion re- 
lation and ours0 the solutions in all of the asymptotic regions (expressions 
[24] and [25] ) are equivalent with the exception of region g and those regions 
for which the transcendental equation of Bogdan et al. [3] does not reduce to 
a simple quadratic (our regions E-H). In order to perform angular integrals 
of the intensity, however, Bogdan et al. [3] make the assumption that the ve- 
locity and material temperature perturbations are independent of angle. This 
does not appear to be a valid assumption in the optically thin limit, since 
coupling between the radiation and material should depend in that case upon 
the direction of the radiation. 



5.3 Lowrie, Morel & Hittinger 

Lowrie et al. [7] analyze perturbations about a nonzero mean flow and focus 
on the initial value problem (solving for u as a function of k) . Under the Ed- 
dington approximation, the dispersion relation ([TBI) is a fifth-order polynomial 
in uj that must be solved numerically. Both for simplicity and due to the fact 
that Kull does not currently currently support periodic boundary conditions, 
we have chosen to focus on the boundary value problem in this work. A pe- 
riodic domain has its own advantages, however, and the initial value problem 
as a numerical test would not be plagued by the transients and reflections we 
have observed when conducting our numerical tests (see $6]). 



6 This is due to the fact that the transcendental equation of Bogdan et al. [3] must 
be expanded to 0(t^ 6 ) due to a cancellation at 0(t®). 
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5.4 Vincenti & Baldwin 



Vincenti & Baldwin [8] perform an analysis similar to that of Bogdan et al. 
[3], with the additional assumptions of negligible radiation pressure and time 
scales much greater than the time scale for coupling between the radiation and 
material (in our notation, r < 1 and r c ^> 1). The correspondence between 
their notation and ours is 

icj «-> T a T^, (27) 

N Bo <-> - (28) 
rc 



and 

N Bu <-> r a . (29) 
We demonstrate in Appendix [B] that their dispersion relation is 



2 

Q6rr c (jt' 2 - r k 2 ) I dfi -/—^ = 0, (30) 



which is equivalent to the transcendental relation 

r- a - r fe - 2 - <16rr c (7r- a - r fc " 2 ) (l - r k tan" 1 r,- 1 ) = 0. (31) 



This reduces to the dispersion relation of Bogdan et al. [3] in the quasistatic 
limit (i.e., the determinant of the matrix at the top of p. 884 of [3 J with s = iuo 
and Cdop = Cdyn = Ctof = 0). Vincenti & Baldwin [8] replace the integral in 
expression (|30|) with a single value of /i = 0.64 and d[i = 0.813 to obtain an 
approximate dispersion relation upon which they base the remainder of their 
analysis: 

(1 - iUrr c ) r fe - 4 + (2.44 - r" 2 + il^rr^ 2 ) r^ 2 - 2.44T a ~ 2 = 0. (32) 



This dispersion relation captures the correct phase speed for the acoustic mode 
in regions a-d and the damping length to within 0.1% (regions a and c), 19% 
(region b) or 20% (region d). For regions A and B, it is the damping length 
that is captured correctly with the phase speeds correct to within 0.1%. For 
regions E and F, the errors in the phase speed and damping length are 19% 
and 27%, respectively. Regions e-g, C-D and G-H are not captured due to the 
additional simplifying assumptions made. 
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5.5 Su & Olson 



The problem considered by Su & Olson [9] is also very similar to the one ana- 
lyzed here. For late times they are essentially the same problem, although the 
analysis in Su Sz Olson [9] is considerably more restrictive in its applicability 
due to additional simplifying assumptions (such as the neglect of hydrody- 
namic motions). After an initial boundary layer in time, their solutions have 
the form of a constant background plus a small perturbation (see their equa- 
tion [31]) The correspondence between their notation and ours is 



I67?", 



167rr c 



(33) 
(34) 



and 



-2 



(35) 



Their equation (21) is a dispersion relation between wavenumber and fre- 
quency. In our notation it is 

(1 - il6-frr c )T^ 2 + 3(1 + ir~ l + lfryr) = 0. (36) 



This captures the radiative diffusion modes in regions B, C, F and G. One can 
show that the velocity perturbation for these modes is usually much smaller 
than the perturbations in temperature and radiation energy, thus validating 
their neglect of hydrodynamic motions!" 8 ] 

Perhaps the most important distinction between our approach and that taken 
by Su & Olson [9] is that their u is imaginary whereas ours is real; as a result, 
their modes do not propagate. With that caveat in mind, their inverse Laplace 
transform operation can be viewed (at late times) as a linear superposition 
of radiative diffusion modes. This superposition introduces two complications: 
1) the high frequency components make numerical evaluation of the semi- 
analytical result difficult, and 2) there is an additional source of error for a code 
comparison due to the attempt to represent a continuum of frequencies with 

7 This implies that their assumption of a cubic temperature dependence for the 
heat capacity is only necessary for the initial temporal boundary layer. 

8 There are portions of region C, particularly for values of e larger than they 
consider, in which hydrodynamic motions become important and this assumption 
breaks down. 
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a discretization. Both of these complications are introduced without testing 
the code beyond what one can do with a single mode. An advantage of the 
approach taken by Su & Olson [5] is that it allows for a wider range of initial 
conditions since it captures the temporal boundary layer. 



6 Numerical Results 



The numerical implementation of the solutions described in the previous sec- 
tion requires driving the boundary of a one-dimensional computational domain 
at frequency u, with material and radiation fluids satisfying the assumptions 
described in £j2j We employ the code Kull for our numerical calculations, a de- 
scription of which can be found in reference [5]. Kull is an Arbitrary Lagrange 
Eulerian (ALE) code, although we only present results with Kull in Lagrangian 
mode. We drive both the radiation temperature and the material velocity at 
one boundary, and we use a Milne boundary condition on the radiation [10] at 
the opposite boundary. Kull does not support an outflow boundary condition 
for the material, so we simply fix the velocity at the opposite boundary and 
stop the calculation before the perturbation reaches the far end of the grid{f] 

The dimensionless measures of temperature and density for a system in LTE 
are a/c and r. For 7 = 5/3 and a mean molecular mass of 0.6, the correspond- 
ing physical temperature and density scales are 

/a\ 2 

T = 4xl0 12 - K (37) 



and 



6 x 10-- Q 



/' 1 gem 3 



3 x 10 14 [ - 

x 



r 1 gem 3 . 



(3f 



We conduct nearly all of our runs with a = 10 4 c, corresponding to a mean 
temperature of 4 x 10 4 K. 

We define our computational domain to be a fixed fraction or multiple of the 
perturbation wave length A. For a given solution to the dispersion relation 
( TT8l) . the opacity is given by 

X = .J^-ii (39) 
XRe[r k ] 



9 Driving at one end with outflow at the other end would be the appropriate bound- 
ary conditions for an Eulerian calculation as well. 



12 



and the driving frequency by 

2vra . , , 

Xr a Re[r k \ 



The numerical length scale can be associated with a physical length scale 
by calculating x based upon a particular frequency-integrated opacity. The 
physical time scale is then determined by this length scale and the speed of 
light. 

Results for the radiative acoustic wave are shown in Figures [3HE1 (regions a-f 
of Figured]). The points are the numerical solutions and the solid lines are the 
analytical solutions. The boundary at which the driving is applied is on the 
left, the computational domain is ten wave lengths, and we run the simulation 
for ten wave periods to ensure that reflection off the right boundary does not 
influence our results. We plot the density perturbation at the end of each 
run. All of these results are at a resolution of 80 zones per wave length. The 
discrepancies at the right hand side of these figures are due to initial transients 
that are not captured by the analytical solution] 10 

Results for the radiative diffusion wave are shown in Figures IMTSI (regions 
A-F of Figure [2]). The computational domain is ten wave lengths in regions 
A, B and D, although we only plot the first wave length since that is the 
length over which the perturbation is damped. We integrate these runs for ten 
wave periods and the number of zones per wave length is again 80. Since the 
damping length in regions E and F is much smaller than a wave length (by a 
factor of ~ 10~ 2 ), we use a computational domain of one half of a wave length 
for these runs, to give ~ 16 zones per damping length. For all of the radiation 
diffusion runs we plot the radiation temperature perturbation at the end of 
each run. 

We include a high resolution result for region C (Figure [TT]) to demonstrate the 
significant numerical cost that can be required to obtain an accurate result. 
The computational domain for this run is one wave length and we integrate 
for two wave periods. The number of zones per wave length is 1600 (20 times 
our nominal value), and the time step was set to the diffusion time scale, as 
opposed to the implicit time integration that was employed for the other runs. 

Figure [13] highlights the fact that coupling to both modes can occur, mak- 
ing it difficult to isolate a single mode. We demonstrate this explicitly in 
Figure [HI which shows Kull results at various resolutions along with the ana- 
lytical solution for both the acoustic mode and the diffusion mode. While we 



The use of outflow boundary conditions would allow these transient features to 
propagate out of the computational domain. 
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are attempting to drive the diffusion mode, it is clear that the acoustic mode 
is being excited as well, with the amplitude of the excitations decreasing as 
we increase the resolution 

Another consideration when comparing the numerical results to the analyt- 
ical solution is that different perturbations dominate in different regions of 
parameter space, and comparisons to the analytical solution are more robust 
when based upon the dominant perturbation. The density and velocity per- 
turbations generally dominate for the acoustic mode, and the radiation tem- 
perature perturbation for the diffusion mode. Another implication of this is 
that care must be taken when setting the overall amplitude of the perturba- 
tions; a small amplitude for the hydrodynamic variables may translate into a 
nonlinear amplitude for the radiation energy density, or vice versa. 

Figures [TpTfT^l demonstrate that flux-limited diffusion does not preserve causal- 
ity for small amplitude waves. Figures [TO] and [T7| are results for an acoustic 
wave (region g of Figure [T]); these were run at a = 0.1c to make this region 
of parameter space more computationally accessible. For a sufficiently small 
amplitude (Figure [TBI) , the numerical solution converges to the analytical solu- 
tion under the diffusion approximation, both with and without a flux limiter. 
A driven wave is superluminal if its wavelength exceeds the wavelength of 
free-streaming radiation at the same frequency; i.e., v p > c for 



The dashed line in Figure [16] shows a wave with a phase velocity equal to c, 
demonstrating the superluminal nature of the excited wave. As the amplitude 
is increased (Figure H7p . the wave begins to steepen into a shock, with the 
flux limiter doing nothing to limit the amplitude of the radiation temperature 
perturbation. 

Figures [18] and [19] show similar results for a diffusion wave (region H of Fig- 
ure [2]). Figure [TS] again demonstrates that the numerical results (with and 
without a flux limiter) are converging to the analytical solution under the dif- 
fusion approximation, which is superluminal in this region of parameter space 
(albeit damped over a wavelength). Figure [19] indicates that as the ampli- 
tude of the wave approaches the nonlinear regime in this case, the flux limiter 
begins to shorten the wavelength. 



The use of outflow boundary conditions on the material would likely reduce this 
effect considerably. 
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7 Summary and Discussion 



We have conducted a systematic comparison of a perturbation analysis of 
the equations of radiation hydrodynamics with the Lagrangian code Kull in 
a wide range of parameter space. We have demonstrated that these solutions 
are a useful benchmark for testing any radiation hydrodynamics code. The 
most important issues to keep in mind when conducting such a test are 1) 
flux-limited diffusion does not capture these solutions in the free-streaming 
limit, 2) coupling to both modes can occur, making comparison with a single 
mode somewhat difficult and 3) care must be taken in setting the perturba- 
tion amplitude, since either the radiation or hydrodynamic perturbations can 
dominate the others by orders of magnitude. 

The primary purpose of a numerical test like the one we have studied is to 
investigate the convergence properties of a numerical discretization. Since our 
focus has been on the test itself rather than on the convergence properties 
of Kull, and due to the complications of transient effects and mode coupling, 
we have not included any convergence results here. We have investigated the 
convergence properties of Kull for most of the solutions shown in Figures I3TTT51 
by calculating the L2 norm of the error between the numerical and analyti- 
cal results (excluding the final wavelength to minimize transient effects). Kull 
converges at second order in most regions of parameter space, as expected, 
particularly when the hydrodynamics and radiation are weakly coupled. In 
regions of parameter space where the coupling between modes is strong, the 
convergence is weaker. A formal convergence test in Kull using the solutions 
described here would require the implementation of outflow boundary condi- 
tions. 

Some additional comments on the breakdown of flux-limited diffusion are in 
order. Figures dS] and [T7| are somewhat of an academic exercise, since the 
radiation energy in this region of parameter space exceeds the rest mass energy 
of the material, and our non-relativistic treatment breaks down [3]. In addition, 
the superluminal propagation we have observed only occurs for waves whose 
amplitude is small compared to the background state, and it is not clear 
that this would have a significant impact on the energy budget of a realistic 
calculation. Finally, we have set up our computational domain to allow the 
excitation of any length and time scales we choose, and even with this freedom 
it was computationally expensive to access the regions of parameter space in 
which the flux limiters break down. A realistic calculation will only be able to 
resolve a small fraction of the length and time scales that we have explored, 
and only with a very large dynamic range will a calculation be able to see the 
excitation of superluminal waves. In any case, when computation has advanced 
to the point where one can easily resolve the length and time scales of free- 
streaming radiation, one should no longer be using flux-limited diffusion. 
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Fig. 1. Parameter space for the radiative acoustic mode under the Eddington ap- 
proximation. The phase speeds are the material sound speed (regions a and b), the 
isothermal sound speed (regions c-e), the radiative sound speed (region f) and c/y/3 
(region g). 
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Fig. 2. Parameter space for the radiative diffusion mode under the Eddington ap- 
proximation. The phase speeds are v p <C a in the white region, a < v p < c in the 
light shaded region and v p = c\/3 in the dark shaded region. 
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Fig. 3. Density perturbation in region a (r = 1CT 3 , r a = 10 4 ) after ten wave periods. 
The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 4. Density perturbation in region b (r = 1(T 5 , r a = 10~ 2 ) after ten wave 
periods. The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 5. Density perturbation in region c (r = 10 3 , r a = 10) after ten wave periods. 
The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 6. Density perturbation in region d (r = 10" 1 , r a = 1CT 2 ) after ten wave 
periods. The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 7. Density perturbation in region e (r = 10 3 , r a = 10 2 ) after ten wave periods. 
The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 8. Density perturbation in region f (r = 10, T a = 10 4 ) after ten wave periods. 
The points are the Kull results, and the solid line is the analytical solution. 
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Fig. 9. Radiation temperature perturbation in region A (r = 10~ 3 , r a = 10 4 ) after 
ten wave periods. The points are the Kull results, and the solid line is the analytical 
solution. Only one-tenth of the computational domain is shown. 
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Fig. 10. Radiation temperature perturbation in region B (r = 10~ 3 , r a = 1) after 
ten wave periods. The points are the Kull results, and the solid line is the analytical 
solution. Only one-tenth of the computational domain is shown. 
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Fig. 11. Radiation temperature perturbation in region C (r = 10~ 3 , r a = 10 4 ) 
after two wave periods. The dotted line is the Kull result, and the solid line is the 
analytical solution. 
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Fig. 12. Radiation temperature perturbation in region D (r = 10 2 , r a = 10 3 ) after 
ten wave periods. The points are the Kull results, and the solid line is the analytical 
solution. Only one-tenth of the computational domain is shown. 
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Fig. 13. Radiation temperature perturbation in region E (r = 3 x 1(T 7 , r a = 3) after 
one wave period. The points are the Kull results, and the solid line is the analytical 
solution. 
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Fig. 14. Radiation temperature perturbation in region E after ten wave periods. 
The solid lines are Kull results with resolution increasing from top to bottom, the 
dotted line is the analytical solution for the diffusion mode, and the dashed line is 
the analytical solution for the acoustic mode. 
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Fig. 15. Radiation temperature perturbation in region F (r = 3xl0" 6 , r a = 3X10" 1 ) 
after one half of a wave period. The points are the Kull results, and the solid line 
is the analytical solution. 
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Fig. 16. Radiation temperature perturbation in region g (r = 10 3 , Ta = l,a = 0.1c) 
after one wave period. The solid lines are the Kull results (both with and without a 
flux limiter; differences are O[10~ n ]), and the dotted line is the analytical solution 
under the diffusion approximation. For reference purposes, the dashed line shows 
a wave at the same driving frequency with a phase velocity equal to the speed of 
light. 
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Fig. 17. Same as Figure [16] with a higher perturbation amplitude. The differences 
between the numerical results with and without a flux limiter are O(10~ 8 ). 
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Fig. 18. Radiation temperature perturbation in region H (r = 10 2 , r a = 1(T 6 ) after 
one wave period. The solid lines are the Kull results (both with and without a flux 
limiter; differences are O[10 -5 ]), and the dotted line is the analytical solution under 
the diffusion approximation. For reference purposes, the dashed line shows a wave 
at the same driving frequency with a phase velocity equal to the speed of light. 
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A Eigenvalues and Eigenvectors 



With d t — > iu; and V — > —ik, the density, pressure and radiative flux pertur- 
bations are given by 



5p t c 5v 

PO T k c' 



(A.l) 



Sp = 8T + Il 5v 
Po T T k c 

and 



where 

A = l + tfr- 1 (A.4) 



and the /'s are flags to keep track of terms. Mihalas & Mihalas [Tf2] set /' = 
everywhere but / = only for the flux in the material momentum equation 
(otherwise / = 1); a consistent treatment has / = /' = 1. For the diffusion 
approximation, / = /' = 0. 

Using the above expressions in the material and radiation energy equations 
gives 

16 -(f-f)-(>- 1 l- 1 T--f) 



and 



-=(l + K 1 + K)—- {A - f'K 1 ) -, (A.6) 
T I c 3A I T V ^ 6 / c ' V ; 



or, equivalently, 
ST 5 



C ~ = iC - l)f - it (A - fir-') * (A.7) 
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where 



.-2 



C=l + ir- 1 + ^. (A. 



Eliminating 5T and 5T r , respectively, from these equations gives 

(lfryr[C7 - 1] + K X C) ^ = ir^ {^A> + [ 7 - ^ (A.9) 



and 

(l6 7 r[C - 1] + ix-'C) f = V f 152 ^^' +T " l) (A.10) 
where 

A' = 1 + (/ - f)vr~\ (A.11) 



The ratio of these expressions, 

ST 16jrA' + 3A( 7 - 1)C 



5T r (167r + ir- 1 )^' + 3A( 7 - 1) 



(A.12) 



demonstrates that the material and radiation temperatures are nearly equal 
in the optically-thick limit (r c ^> 1 and ^> 1). 



The material momentum equation is 

3i( 7 -l)y c ' T ' 3A(j-l) T ( 



™ (V - - ^ + £ + r^TT^ = 0, (A.13) 



where 

A^l+i/r" 1 (A. 14) 



has been defined to indicate that Mihalas & Mihalas [Tf2] ignore the frequency 
dependent term in this equation (for them, / = 0). Replacing 5v with expres- 
sion (IA.9I) and 5T r with expression flA.12[) gives the dispersion relation: 

c 4 T- fc " 4 + c 2 r fc ~ 2 + c = 0, (A.15) 
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with 

C4 = 1 — H6rr c , 



(A.16) 



c 2 = 3A(l + it- 1 ) - t- 2 (1 - il6 7 rr c ) 



and 

Co = -3Ar- (1 + 16 7 r + ^r 1 ) (l + f ^"[y ) ■ ( A ^) 

Under the assumptions of Mihalas & Mihalas [Tf2] (/' = / = and / = 1), 
these become 

c 2 = 3 (l + ir" 1 ) 2 - r" 2 (1 - zl6 7 rr c ) 

and 

c = -3 (l + ir' 1 ) r~ 2 (l + 16 7 r + ir' 1 ) , (A.20) 

with C4 unchanged. These coefficients match those of expression (3.12) in Mi- 
halas & Mihalas [lj. The eigenvector relationships are 

_ 1 \ 5T . _j /16 7 r , \ 5f 



and 



- 1] + «r-'C) £ = ir- + [7 " 1]C) ^ (A.21 



(IMC - 1] + ir-'C) ^ = ir t - feJSl + 7 - l) ( A.22 



A self-consistent treatment under the Eddington approximation has f' = f 
f — 1, which gives 



c 2 = 3 + ir c l ) - r a 2 (1 - 2i6 7 rr ( 

6 7 r + zr c 
3(7-1) 



+16rf5 + 3ir B - 1 + 1 ^ r + i :n ( A - 2 3) 
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and 



c = -3r- 2 (l + 16 7 r + tr' 1 ) (l + ir" 1 + ^^y) ■ ( A - 2 4) 



The self-consistent eigenvector relationships are 

(lMC _ 11+jr - lc) | = fcsrl (_^_ + [7 _ 1]c )^ (A , 5) 

and 

The diffusion approximation (/' = / = / = 0) gives 



c 2 = 3 fl + «t c M - r a 2 (f - HQ^yrr, 

6 7 r + ir c 
3(7-1) 



+ 16r f 5+ 16 y + ^^ (A.27) 



and 

c = -3r" 2 (l + 16 7 r + ir" 1 ) . (A.28) 

The eigenvector relationships under the diffusion approximation are given by 
expressions (IA.21[) and (IA.22I) with C replaced by 

G 1 = 1 + zr- 1 + (A.29) 

o 



We have derived the approximate expressions n2M and ( 1251) both analytically 
and by a semi-empirical approach described below. The positive and negative 
branches of the dispersion relation are given approximately by 

^*-f(l + f) (A.30) 

and 
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where the second term in parentheses is only required to derive the expressions 
in regions a and c. Asymptotic expansions are performed on the coefficients 
of the dispersion relation and the leading terms are inserted in the above 
expressions. As it turns out, these asymptotic expansions are not always easy 
to perform, and we have been unable to derive the approximate expression 
in region g analyticallyP^I A more straightforward semi-empirical approach is 
to calculate linear fits of the full solutions on a log-log plot to determine the 
scaling of the solutions with the various parameters. We have done this to 
obtain the expression for region g as well as to check our analytical approach 
for the other regions. 



B Details on the Vincenti & Baldwin Analysis 



Equation (51) of Vincenti & Baldwin [8J is 



m] + H " [ ^ ]) = ( B ~ m " rlH " [0] ) E2{Nbu ^ 

' E 2 {N Bu \i - £]) (#'[£] + r'H'i'ii}) di 



oo 



+ / E 2 (N Bu [i - £]) (#'[£] + ^H"^}) dl (B.1) 



where if is a dimensionless perturbation amplitude, B is a boundary condi- 
tion, a prime denotes a derivative with respect to £ (a dimensionless spatial 
coordinate) and 

i 

E 2 {z) = f e z/s ds. (B.2) 



At this point Vincenti & Baldwin [8] replace H with HjCje c ^ (a sum of expo- 
nentials) and approximate E 2 with a single exponential. This approximation 
is unnecessary, however. Making only the former substitution gives 



N, 



Bo 



Bu 



(l + cJ)C 



B 



) E 2 (N Bu O 



12 An important consideration with the analytical approach is to only perform 
asymptotic expansions of complex expressions that appear in the numerator; com- 
plex expressions in the denominator should be converted to real expressions via 
their complex conjugate. 
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£ i 



-Jdijdn e-^S/v+^^ZjCj (l + 7- 1 c|) Cje c ^ 
o o 

oo 1 

+ Jd^Jdfi e NB »V ft - NB »& ,i E j c j (l + 7 _1 cf) Cje **. (B.3) 



£ o 



Changing the order of integration, taking the £ factor out of the integral over 
£ as well as the constants out of both integrals, and reversing the limits of 
integration on the second integral over £ gives 



N Bo 



Bu 



Ej (l + cj) Cje^t = (B- T, 3 C 3 [l + t- 1 ^]) E 2 (N Bu O 



-Z jCj (l + 7-^) Cj 



( 1 i 

J dfi e- NB «V» J di e ( c i +J WM)£ 

\o 
1 £ 
+ J dfi e NB »* /ft J di ete-Wd* 
00 



Performing the integrals over £ gives 



(B.4) 



iV, 



Bo 



87^Bn 

1 

-E jCj (l + 7~ 1 c 2 ) J dfi 



Ej (l + cfj = {B- EjCj [1 + 7 - 1 c J 2 J) E 2 (N Bu O 

cjfi + N Bu cjfi - N Bu 



1 

-EjCj (l + 7~ 1 c 2 ) Cj J dfi 



Cjfi + N Bu 



(B.5) 



where Re(cj) < has been assumed (this is necessary for the perturbations 
to remain finite as £ — > 00). The integral over /1 in the second line above is 



dfi 



2 2 



C; 



1-^tanlW °> 



N Bu 



(B.6) 



The terms proportional to e c ^ in expression flB.51) yield the dispersion rela- 
tion^ 



, , 2 . lQN Bu 
1 + c, - z— — 



(7 + c?) 



tanh 1 



iV, 



Bu 



(B.7) 



13 This is equivalent to expression (|31|) in the text after making the substitutions 
OH)- (1291). 
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while the remainder of the terms imply 



± 

fdf* 



N. 



Bu 



cjn + N Bu _ 



(B.8) 



This expresses the boundary condition at the wall. It can be rewritten in 
terms of exponential integrals and integrated, but it does not appear that this 
expression can be satisfied for all £. This may be due to the fact that the 
discrete spectrum is not sufficient to match the particular boundary condition 
chosen by Vincenti & Baldwin [8]. Equation flB.71) also gives rise to a contin- 
uous spectrum of modes when the wavenumber is purely imaginary and the 



optical depth is unity (cj/N Bu = ±1)1 I Including the continuous spectrum in 



the decomposition of H(£) would likely alleviate the discrepancy in equation 

id. 
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This corresponds to a singularity in the dispersion relation when the argument 
of tanh -1 is ±1. See Bogdan et al. [3] for further discussion. 
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